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Abstract 

Let Pn be the probability for a planar Poisson- Voronoi cell to have 
exactly n sides. We construct the asymptotic expansion of log pn up to 
terms that vanish as n — > oo. Along with it comes a nearly complete 
understanding of the structure of the large cell. We show that two in- 
dependent biased random walks executed by the polar angle determine 
the trajectory of the cell perimeter. We find the limit distributions of 
(i) the angle between two successive vertex vectors, and (ii) the one 
between two successive perimeter segments. We obtain the probabil- 
ity law for the perimeter's long wavelength deviations from circularity. 
We prove Lewis' law and show that it has coefficient 1/4. 

1. Introduction 

Cellular structures are observed in a great diversity of situations in science 
and technology. A common phenomenological description of such structures 
makes use of Voronoi cells PP , constructed around centers which for simplicity 
are referred to as "particles": the Voronoi cell of a particle a is made up of 
those points in space for which a is the closest particle. Hence Voronoi 
cells constitute a "tessellation", i.e. a division of space into nonoverlapping 
regions. An excellent review of tessellations, with in its introduction many 
historical and anecdotal details, is due to Okabe et al. It offers a thorough 
discussion of much fundamental and applied work in this area. This includes 
applications to biological tissues, mesh generation in numerical algorithms, 
coding in telecommunication, etc. 

Whereas the Wigner-Seitz cell associated with the regular lattices of solid 
state physics is a special case, the generic Voronoi tessellation is based on 
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randomly distributed particles; for uniformly distributed point-particles the 
result is called a Poisson- Voronoi tessellation. 

Voronoi cells appear in physics for many different reasons. Meijering |3] 
used them to characterize the end result of a crystal growth process starting 
from independent centers. In the theory of the solid-liquid transition they 
serve to define nearest-neighbor relations between atoms, which in turn allow 
the identification of lattice defects. Recently Voronoi tessellations of fractal 
sets have been considered |3]. In field theory the use of random lattices and 
the associated Voronoi tessellation Ej was motivated by their advantage of 
being statistically invariant under arbitrary translations and rotations while 
still preserving a short-distance cutoff. 

Here we are interested in two-dimensional Poisson- Voronoi tessellations. 
Planar Voronoi cells are convex polygons. The first and foremost quantity 
that appears in both theoretical and applied work is the probability p„ for 
a planar Voronoi cell to be n-sided. Approximate values for have been 
tabulated by many authors j^I; they peak at n = 6 and fall off rapidly 
for larger n. Connected questions concern the statistics of the cell area, its 
perimeter, and its angles, as well as correlations with neighboring cells. They 
have been studied by methods ranging from heuristic arguments and Monte 
Carlo simulations to exact methods. Whereas some statistical properties 
are readily obtained, others, among which the fraction p„, still defy solution 

0121. 

Interest in planar Voronoi cells having a large number n of sides first 
arose in biology. Lewis' empirical law [71 12] , formulated more than seventy 
years ago, states that when n becomes large, the average area of the n- 
sided cell grows as ~ c(n — 2)/z/, where v is the particle density. In spite of 
attempts, this asymptotic proportionality with n has not so far been proved 
[H]. Numerical estimates of the coefficient c are in the range 0.199 — 0.257 
[HI El 121 • In metallography large cells were observed by Aboav ^1] in the 
arrangement of grains in a polycrystal. His empirical law states that the 
neighbor of an n-sided cell has itself on average + ain~^ neighbors. Only 
recently Hug et al. jT2| showed rigorously that when a Voronoi cell becomes 
large (in an appropriate sense), its shape tends towards a circle. In a different 
approach mathematicians study those large cells that allow for an inscribed 
circle whose given radius becomes asymptotically large jT^]- 

Drouffe and Itzykson (DI) estimated the probability Pn for a planar 
Voronoi cell to be n-sided by Monte Carlo simulations for n up to 50. The 
rejection rate limits the statistical precision of generating such very rare 
events. All other simulations cited in Ref. (see e.g. pUj) were restricted 
to n ^ 14, so that the work by DI is still today the main reference for the 
large n behavior. 

This note is the first announcement of a collection of new analytical results 
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on planar Poisson-Voronoi tessellations. All of them are by-products of our 
asymptotic evaluation of p„, and together they lead to a nearly complete 
understanding of the large n-sided cell. Our analysis uses an accumulation of 
classical mathematical methods; however, the asymptotics is beset by a great 
number of difficulties, of which the correct choice of variables of integration 
is merely the first one. It will be possible here to indicate only the main steps 
and the main results. 

We begin by a precise statement of the question. Let a subdomain of 
the plane, of area L^, contain N particles with tags a, placed at positions 
fHa chosen independently and with uniform probability, and let N, L ^ oo 
at fixed N/L'^ = v. We consider the Voronoi cell of an arbitrarily chosen 
particle ag whose position will be the origin of the coordinate system. The 
perpendicular bisectors of the fH^ (with a ^ ao) pass through the midpoints 
Ra = |9^a5 which are uniformly distributed with density Av. Each side of the 
cell is a segment of one of the perpendicular bisectors. An ra-sided cell may be 
fully specified by n midpoint vectors Rm (m = 1, . . . , ra); or alternatively by 
the n vertex vectors Sm {m = 1, . . . , n) that connect the origin to its vertices. 

The probability Pn has been expressed jHl El HI] as the 2n-dimensional 
integral 



where we scaled lengths such that 4z/ = 1; it is analogous to a partition func- 
tion in statistical mechanics. The functions x and A enforce two conditions 
necessary for the n vectors Rm to define a valid Voronoi cell. First, each of 
the n associated bisectors must contribute a nonzero segment to the perime- 
ter. Calka formulates this condition as a set of n inequalities, and x is 
the indicator function of the phase space domain where they are satisfied. 
Secondly, the union of the n disks centered at and of radius must 
be void of midpoints. The probability for this is e""^, where A is the area 
of the void region. Explicit expressions for A were given by DI [Hj and by 
Calka ^3]. Eq. has been rewritten and evaluated to seven decimal places 
for ra = 3, but has otherwise remained untractable analytically. 

2. The perimeter as a Markov process 

Our approach originated from an attempt to describe the perimeter of 
a Voronoi cell as a Markov process with the polar angle as independent 
variable. We order the midpoints Rm = {Rm, ^m) such that < $i < 
. . . < $„_i < $„ = 27r and denote by Sm = {Sm, ^m) the vertex where the 
bisectors through Rm-i and Rm intersect. Let Pm be the angle between Rm-i 
and Sm, and jm the one between Sm and Rm- These angles may be positive 
or negative, since a midpoint need not (and for large n generically will not!) 
itself be part of the side that it defines. The Markov process description 




m=l 



(1) 
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leads to a transition probability between two successive midpoints which has, 
essentially, the form Rl^Tme'-^"". Here Am, which satisfies A = Ylm=i-^"^^ 
is the area of a region delimited by three circular arcs that should be void 
of midpoints if the perimeter is to travel uninterruptedly between Sm-i and 
Sm] and Tm = cos7mSin(/3m + 7m)/cos^/3m is the (unnormalized) probability 
density for the perimeter to turn by + 7^ at vertex S^- Considerable 
algebra then gives 

Pn = ^ I dpdj6{ J2{Pm + lm)-2n) 6{G) 

m=l 

r " n r°° 

x\Up"mTm] / di?i?^"-V-«^^(^'^), (2) 
which calls for several explanations, (i) We abbreviated 

-tt/2 pTT/2 p-7T/2 ptt/2 



d/3d7 = / dp^ d7i / dP2 / d72 . . . 

J-7r/2 J-fS^ J-li ''-P2 

/tt/2 pTT/2 
d(3n d7„^(7„ + /3i). (3) 

Because of how they are nested, these integrations impose an orientation 
on the perimeter, (ii) For an arbitrary set of angles (/?, 7) = {Pm,1m} the 
perimeter will spiral instead of close onto itself after a full turn of 27r; the 
factor S{G), where G = (27r)~^ Ylm=i log(cos7m/ cos /5m), enforces its closure, 
(iii) The ratios pm = Rm/ R, where R is the average midpoint distance R = 
YlZi=i functions of the angles defined by 

cos7mCOS7m-i . . •COS71 , . 

Pm = n n IT Pn [V 

COS Prn COS Pm-l ■ ■ ■ COS Pi 

for m = 1, . . . , n — 1, together with the sum rule Ylm=i = 1- (i^) 
Finally, we set A = tcR^A{P,'j). 

Eq. (j2|) fully defines our starting point. As a check one may derive it 
directly from by a coordinate transformation, in which the factors p^Tm 
then appear as the Jacobian. The decisive advantage of Q over ((H) is that 
the integration limits have been made explicit so that the function x is no 
longer needed. 

Integrating on i? in yields |(n — 1)! [7ryl(/5, 7)]"". If, as will appear, 
A{i3,'~f) differs negligibly from unity for n —>■ 00, a steepest descent analysis 
of the same integral shows that it draws its main contribution from R within 
a width of order n° around the saddle point Rc = ^/n/n. 

Let now Pn(^) be expression Q but with the 27r in the delta func- 
tion replaced by a new variable Q, and let Pn{s) be its Laplace trans- 
form with Laplace variable ns. Here we pass to a "canonical ensemble" 
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of angles with weight exp[— ns X]m=i(/^m + 7m)]- We define EI(s) by writing 
Pn{s) = J dpdj e"^'-*''. The evaluation of this integral constitutes the real 
problem: it is on all possible shapes {P, 7) of a Voronoi cell with average 
midpoint distance scaled to unity. 

3. Asymptotic expansion 

Anticipating that we will be able to identify a zeroth order problem we 
write H(s) = Ho(s) + V. Hence 

d/?d7 e-i^^^) = (e-^)„ / d/3d7 e-=°('') , (5) 



where (. . .)(, is the average with weight exp[— E[o(s)]. The appropriate choice 
of Ho, if there is one at all, depends delicately on how the variables of in- 
tegration in (jSj) are assumed to scale with n. We certainly expect that for 
large n only small Pm and 7^ will contribute, but the exact scaling of these 
angles with n is not a priori evident. The calculation bears out that 



71" " JJ (/5m + 7m) exp[-ns{Prn + 7m)] (6) 



is the right choice for Hq. 

We now transform to another set of angles, viz. C,m = Pm + 7m and 
^?m = 7m + Pm+1 whcre m = 1,. . . ,n and /5„+i = Pi. Clearly C.m is the 
angle between two successive midpoint vectors and r]m the one between two 
successive vertex vectors. Since the set {^m, Vm} does not fix the relative angle 
between the systems of midpoint and of vertex vectors, we complete it by 
one of the original angles. Pi. Setting = Yli=i ^m = Pi + J2T=i Ve 

we obtain for m = 1, . . . , n the inverse relation 

/9m = ^m - $m-l, 7m = $m " ^m • (7) 

The "zeroth order" problem is constituted by the integral that is left in 
the RHS of (0) when one sets V = 0. It can be solved when reformulated 
in terms of the and r]^. Since ^m,Vm > and J2m^rn = J2mVm = 2n, 
they must scale as = n~^Xm and rjm = n~^ym- In the limit n ^ oo 
the Xm and Um niay be integrated from to cxd and the problem of the 
nested integrations disappears. One finds /d/5d7 exp[— EIo(s)] = (3n — 
2)! 7r^"(ns)~^"^^/(2r?,)!, up to contributions that vanish exponentially with 
n. Let Pn\^) be the Laplace inverse of this zeroth order expression; one 
readily finds Pi°''(27r) = (87r^)"/[47r^(2r2)!]. The inverse Laplace integral has 
a saddle point at s = Sc = 3/^2; if we can show that for n — > cxo the prefactor 
(exp(— V))o in Eq. © has a finite nonzero value in s = Sc, then it follows 
that pn = (exp(-V))3/2.i'^°^(27r). 
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Figure 1: Schematic picture of the two random walks associated with the perime- 
ter of a Voronoi cell. The walk joining tliG vertices S'lfi is the actual perimeter; the 
other walk joins the midpoints Rm- In the limit n ^ oo the increments of the polar 
angle along each walk are independent random variables distributed according to 
Eq. The perimeter segments are on the perpendicular bisectors (dashed lines). 

We pause for a few comments. The zeroth order calculation has a useful 
corollary: in the limit n — >• oo the scaled angles Xm and i/m become indepen- 
dent random variables with laws u{xm) and v{ym) given by 

u{x) = 7^^\e-"/^ v{y) = {2n)-'e-y/^^^^ . (8) 

It follows, again for n ^ oo, that {$rrt}m=o {^m}m=o ^^'^ mutu- 
ally independent one-dimensional random walks, each with a drift of 27r/n 
per step and with independent increments; they describe how the polar co- 
ordinate of the midpoint vector Rm = {Rm, ^m) and of the vertex vector 
Sm = {Sm,'^m), respectively, increases as one progresses along the perime- 
ter. These angular random walks drive the radial coordinates Rm and Sm and 
couple them together, which leads to the double-stranded two-dimensional 
walk depicted in Fig.[T] 

Let us now pursue the analysis. The validity of the zeroth order results 
hinges on our ability to prove that the first factor in Eq. (jHJ, (exp(— V))q, 
exists and is subdominant in the large n expansion of logp„. Actually, ev- 
erything has been so prearranged that this factor is of order Moreover, 
it has an interesting structure. 

The regular n-sided polygon is an obvious point of symmetry in phase 
space. DI obtained an exact lower bound for pn by expanding the vertex 
vectors Sm about their regular polygon positions. This procedure preserves. 
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however, the long range polygonal order. Here we will advance as in the 
theory of elasticity, where one considers the deviations not of the positions, 
but of the interatomic distances, from their ground state values. Hence we 
expand V in powers of 5xm = Xm — x and 5ym = Vm — y (where x = y = 2n). 
This procedure allows, in principle, for large deviations from the regular 
n-sided polygon. 

Terms from various sources contribute to V and have to be examined one 
by one. The scaling of the and rj^ and their asymptotic independence 
determine how the other variables scale; Eq. ((Zj) implies, in particular, that 
the original angles scale as Pm = n~2b^ and 7^ = n~2Cm,- One finds that 
V = Vo + n-5Vi + . . . . Hence (exp(-V))o = (exp(-Vo))o [1 + 0{n-^)]. It 
appears that Vq is quadratic in the 6xm and 6ym with contributions coming 
from (i) the 0{n~^) terms in the expansion A = 1 + n~^Ai + . . ., and (ii) 
the product Y[m=i Pm'^rn/ {Pm + 7m)- In tcrms of the Fourier components 
Zg = (27m2 )~i J2m=i ^'^^^^^^^^m-, where z = x^y and q is integer, one finally 
finds, with s = 3/(27r), 

^o = J2 [(^'' + 2r')(£. - yg)i^-g - y-q) + 2g-'(x, - y,)y.,] , (9) 

where q runs through all nonzero integers. The decay of the Fourier coeffi- 
cients as ~ q~^ for large q is characteristic of interacting Coulomb charges. 
Since the Xq and yq are sums of independent variables, they are Gaussian dis- 
tributed with variances determined by (jH))- We can show that the quantity 
exp(— Vo) may be averaged by Gaussian integration: in terms of Coulomb 
charges this amounts to a Debye-Hiickel approximation, which becomes exact 
in the high temperature limit; here the inverse "temperature" is 1/n and the 
exactness of our procedure follows from a Hubbard-Stratonovich trans- 
formation. Gaussian integration then leads to (exp(—Vo)) 3/2,1 = C where 
C = (-'- ~ + 4(3'~^)~^ = 0.344... Hence the finite size correction fac- 
tor C is the partition function of the deformational modes of the perimeter 
with wavelength 27r/g, the main contributions coming from the long wave- 
lengths. 

4. Results 

Upon combining the zeroth and first order result we conclude that the 
probability p„ for a Voronoi cell to be n-sided behaves as 

which may also be written as the asymptotic equality 

logp„ - -2nlogn + nlog(27rV) - ^\og{2^7r^C-^n). (11) 
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Figure 2: Straight solid line: the asymptotic formula, Eg. (|TT|) . Open squares: 
Monte Carlo data by Drouffe and Itzykson (Ref. |H], Table 1). Dashed line: 
Eq. (Illj) with an additional term an~2, where we adjusted a = —0.24 to 
fit the data. 

We compare this to existing work. DI assume an expansion similar to (jllj] 
and provide strong analytic an numerical evidence for the coefficient of the 
nlogn term to be equal to —2. This has now been confirmed. Further- 
more, DI prove that the term proportional to n has a coefficient larger than 
log(7r^/e) = 1.289..., and from their numerical data estimate that it is ac- 
tually between 4 and 5. Our value log(27r^e^) = 4.982... satisfies DTs exact 
inequaUty, as of course it had to, and is just inside the range of values that 
these authors judged likely. 

In Fig. 121 our asymptotic result Eq. (fTT|) is shown as a straight solid line of 
slope — |. It is compared to the Monte Carlo data of Ref. |8j (open squares) 
in the range 50 > n > 7 (the abscissa value of 2 corresponds very nearly to 
n = 7). The first term not displayed in the asymptotic series of Eq. (jlip will 
be of the form an~2; if, venturing at this point beyond our analytic results, 
we choose its amplitude to fit the data, very good agreement (the dashed 
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line) with the Monte Carlo results is obtained for a = —0.24 ± 0.02. We 
believe that the downward trend in the DI data for n ^ 30 (abscissa ^ 0.75) 
is due to the great difficulty of such simulations; in that region the authors 
estimated their error bars to be at least of the order of the data themselves. 

Consistency of the asymptotic scaling with Eq. (@]) requires that = 
1 + n~2rm, where remains of order n° as n — > cx). Hence the midpoint 
distances Rm = Rpm deviate by order n° from their average i?, which itself 
is typically within nP from R^. Therefore the area of the n-sided Voronoi 
cell is sharply peaked around an average equal to vri?^ = n/(4z/), where we 
restored the particle density. This demonstrates that to leading order Lewis' 
law holds. Its coefficient c = \ agrees with DI's observation that c ~ |. 

It is amazing that the analysis of a problem defined on only a random 
set of points in the plane evokes associations with so many phenomena in 
statistical physics. A full account TB^ of this work is in preparation. In an 
extension of it we will revisit Aboav's law. 
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